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ABSTRACT 

An  existing  computer  program  for  the  geometrically  non- 
linear analysis  of  arbitrarily  loaded  shells  of  revolution 
was  modified  to  incorporate  initial  geometric  imperfections 
into  the  buckling  solution  of  an  axially  loaded  circular 
cylindrical  shell.   Several  different  boundary  conditions 
were  considered,  and  the  predicted  buckling  loads  were 
examined  to  determine  the  significance  of  the  initial  imper- 
fections and  the  boundary  conditions  for  the  shell  stability 
problem.   The  computed  buckling  loads  were  compared  to  those 
found  experimentally  by  Arbocz  and  Babcock  at  the  California 


generally  higher  than  the  experimental  results.   It  was 
indicated  that  imperfection  sensitivity  is  dependent  on  the 
boundary  conditions  of  the  cylinder. 
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I.   INTRODUCTION 

The  study  of  the  instability  of  thin,  isotropic,  circular 
cylindrical  shells  has  long  been  troubled  with  a  lack  of 
correlation  between  the  theoretically  predicted  buckling 
loads  and  the  experimental  results;  the  experimental  values 
of  buckling  are  generally  much  lower  than  the  theoretical 
solutions.   This  lack  of  agreement  has  been  attributed  to 
imperfection  sensitivity,  i.e.,  inherent  imperfections 
in  the  shape  of  the  shell   which  are  unaccounted  for  in  the 
classical  buckling  theory   trigger  the  buckling  modes  pre- 
maturely, and  to  the  use  of  incorrect  boundary  conditions 
in  t:he  analysis,  _l  .  e  ,  Llie  b^undcii."^  conditions  used  in  the 
analysis  are  not  equivalent  to  those  of  the  test  specimen  [1], 

The  recent  development  of  several  sophisticated  digital 
computer  programs  for  the  geometrically  nonlinear  analysis  of 
shell  structures  has  made  it  possible  to  include  both  the 
initial  imperfections  and  arbitrary  boundary  conditions  in 
the  analysis.   For  example,  a  digital  computer  program  for 
the  geometrical  nonlinear  analysis  of  arbitrarily  loaded 
shells  of  revolution  has  been  developed  by  Ball  [2].   A 
modification  has  recently  been  made  to  this  program  to 
include  the  initial  imperfections  in  the  analysis.   However, 
this  modification  does  not  eliminate  all  of  the  previous 
analytical  difficulties  because  in  order  to  predict  the 
buckling  load  of  a  particular  shell,  the  imperfect  shape  of 


that  shell  must  be  known,  and  the  actual  boundary  conditions 
of  the  specimen  must  be  known.   With  regard  to  the  first 
problem,  Arbocz  and  Babcock  [3]  have  developed  a  means  of 
mapping  the  topography  of  cylindrical  shells  and  can  thus 
describe  the  deviations  from  a  perfect  circular  cylinder. 
The  objective  of  this  thesis  is  to  use  the  imperfection  data 
for  one  of  the  cylinders  presented  in  Ref.  [3]  and  Ball's 
modified  computer  program  to  obtain  stability  solutions  for 
several  boundary  conditions.   The  buckling  loads  are  compared 
with  the  experimental  results  presented  in  Ref.  [3]  and  are 
compared  with  each  other  to  determine  the  effect  of  the 
various  boundary  conditions. 


II.   DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

The  computer  program  used  in  this  study  is  capable  of 
analyzing  thin  isotropic  shells  under  the  following  condi- 
tions : 

1)  The  geometric  and  material  properties  must  be  axi- 
symmetric,  but  may  vary  along  a  meridian. 

2)  The  initial  imperfections,  pressures,  and  tempera- 
tures must  be  symmetric  about  a  meridional  plane. 

3)  The  boundaries  may  be  free,  fixed,  or  elastically 
restrained. 

The  governing  partial  differential  equations  are  based 
on  Sanders'  nonlinear  shell  theory  for  the  condition  of 
small  strains  and  moderately  small  rotations  (Appendix  A) . 
Inplane  and  normal  forces  are  accounted  for  but  rotary 
inertial  forces  are  neglected.   The  set  of  partial  differ- 
ential equations  is  reduced  to  an  infinite  number  of  sets 
of  four  second-order  differential  equations  in  the  meridional 
and  time  coordinates  by  expansion  of  the  dependent  variables 
in  a  sine  or  cosine  series  in  the  argument  n9  where  6  is 
the  circumferential  coordinate  and  n  is  the  Fourier  index. 
Nonlinear  coupling  terms  are  treated  as  pseudo  loads  and 
trigonometric  identities  are  used  to  uncouple  the  sets  of 
equations.   Derivatives  with  respect  to  the  meridional 
coordinates  are  replaced  by  the  conventional  finite  differ- 
ence approximation.   This  leads  to  sets  of  algebraic 


equations  in  the  four  dependent  variables  U    ,  V    ,  W    , 
and  M    where  U,  V,  and  W  are  the  displacements  in  the 
meridional,  circumferential,  and  normal  directions  respec- 
tively, and  M   is  the  meridional  bending  moment  per  unit 
length.   At  each  load  or  time  step,  an  estimate  of  the 
solution  is  obtained  by  extropolaticn  from  solutions  at 
the  previous  steps .   The  sets  of  equations  are  repeatedly 
solved  using  Gaussian  elimination,  and  the  pseudo  loads 
are  recomputed  until  the  solution  converges.   If  convergence 
is  not  achieved  in  a  specified  number  of  iterations,  the 
load  is  reduced  by  a  factor  of  five.   If,  after  a  specified 
number  of  load  reductions,  the  solution  fails  to  converge 
the  problem  is  terminated.   If  the  load-displacement 

b'         --'  -  -C       •!..        nUal  1        -—       —  c      #-V»  a       --f  J —-."«,»      i_.,v-.^  knnlrl  inn 

CUOvlui       Ox.        UlC       oucj-a        a-^>       wj-        ^jli^      _>^j-»_^.n.-i.±.';j        ._  j  j-^  ,       /Juc-J\j.i.»  j 

is  presumed  to  have  occurred  at  the  final  load. 

The  original  program  described  in  [2]  was  not  designed 
to  solve  shell  stability  problems  where  initial  imperfections 
had  to  be  accounted  for.   However,  it  soon  became  apparent 
that  such  a  capability  was  very  desirable  and  the  program 
was  modified  accordingly.   The  equations  and  FORTRAN  state- 
ments that  include  the  imperfection  terms  are  presented  in 
Appendix  A. 
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III.   BOUNDARY  CONDITIONS 


In  the  analysis ,  the  boundary  conditions  for  the  cylin- 
der are  on  N   or  U,  N  fl  or  V,  Q   or  W,  and  8W/8s  or  M  , 
where  N   and  N  fi  are  membrane  forces  per  unit  length,  and 
Q   is  the  transverse  force  per  unit  length.   The  boundary 
conditions  are  represented  by  the  matrix  equation 


f  "\ 

N 
s 


[«]  K 


N 


aw 

8s 


r     "\ 

U 


>+  [A]  K 


V 
W 

M 

s 

V.  J 


\    =     U> 


where  [ft]  and  [A]  are  4x4  matrices  and  {&}  is  a  1  x  4 

column  matrix.   There  is  one  such  equation  for  each  boundary 
of  the  shell.   Note  that  the  analysis  assumes  that  the 
boundary  conditions  are  the  same  for  all  values  of  n ,  where 
n  is  the  mode  number  in  the  circumferential  direction. 
The  testing  configuration  used  in  Ref.  [3]  for  the 
axially  loaded  shell  is  shown  in  Fig.  1.   The  ends  of  the 
shell  are  attached  to  an  end  ring  and  a  load  cell  with 
Cerrelow,  a  low  melting  point  alloy.   The  ends  of  the  shell 
are  clamped  such  that 


W  =  3W/3s  =  V  =  0 


(1) 


The  axial  load  on  the  test  specimen  was  applied  by  means 
of  four  screws  located  at  the  corners  of  the  square  end 
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End  Rate- -4»-     /  j^C(/^//t//\ 

Devcon 
End  Ring 
Test  Specimen 


Cerrolow 


Load  Coll 

Devcon 
End  Plate^ 


FIG. I     CYLINDRICAL    SHELL  TESTING    CONFIGURATION 
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plates.  At  each  load  increment,  the  screws  were  adjusted 
to  give  as  uniform  a  circumferential  load  distribution  as 
possible . 

In  order  to  investigate  the  problem  completely ,  several 
boundary  conditions  were  formulated  for  the  computer  program, 
These  formulations  follow  below. 

If  the  axial  load  is  axisymmetric,  then  on  each  boundary 


and 


N(0)  -H 
s      o 


U(n)  *  0 


N(n)  =  0 
s 


n  =  1,2,.. 


(2) 


n  =  0,1,2, 


This  condition  will  be  referred  to  as  boundary  condition  A. 
Note  in  Eq .  (2)  that  the  condition  on  N   is  a  function  of  n 
This  necessitates  changes  in  the  program  to  account  for 
different  values  of  the  first  element  of  I    for  n  =  0  and 
n  j-    0.   These  changes  are  given  in  Appendix  B. 
Another  possibility  is 


N(0>  -N 
s      o 


N(n)  =  0 
s 


n  =  1,2,... 


at  one  end  of  the  cylinder,  and 


u(n)  =  0 


,    N^V  0 
s 


n  =  0,1,2,... 

at  the  other  end.   This  condition  will  be  referred  to  as 
boundary  condition  B. 

A  third  interpretation  of  the  test  set-up  consists  of 
rigid  end  plates  incrementally  shortening  the  cylinder. 
This  is  represented  at  the  boundaries  by 
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U(0)  =  U     ,    U(n)  =0        n  =  1,2,... 
o 

at  one  end,  and 

U(n)  =0  n  =  0,1,2,. . . 

at  the  other  end.   This  condition  will  be  referred  to  as 
boundary  condition  C. 

For  purposes  of  comparison,  a  simply  supported  cylinder 
will  also  be  analyzed  using  boundary  conditions  A,  B,  and 
C.   For  this  case,  equation  (1)  becomes 

W  =  M  =  V  =  0  (3) 

s 
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IV.   INITIAL  IMPERFECTIONS 

Arbocz  and  Babcock's  topography  mapping  procedure  is 
described  in  detail  in  Ref.  [3].   The  final  imperfection 
data  consisted  of  deflections  from  a  perfectly  circular 
cylinder,  denoted  here  by  W,  around  fifteen  circumferences 
spaced  0.5  in.  apart.   The  first  and  the  last  circumferen- 
tial mappings  occurred  at  a  distance  of  0.125  in.  from  the 
ends  of  the  shell.   Forty-nine  data  points  were  taken  around 
each  circumference  with  the  first  and  the  last  being  the 
same  point  on  one  arbitrarily  chosen  meridian. 

As  noted  previously,  Ball's  program  requires  a  meridional 
plane  aoour  which  the  imperfect  ^liuJci  ii  symmetric.   Thic 
requirement  is  not  generally  met  by  imperfect  cylinders; 
thus  symmetry  must  be  created  in  some  optimum  way.   An 
investigation  to  determine  the  degree  of  symmetry  was  con- 
ducted by  generating  a  cosine  Fourier  series  for  each 
circumference  using  every  meridian  of  data  as  a  base 
meridian.   In  all  cases,  the  phase  angles  for  each  circum- 
ference varied  considerably,  thus  eliminating  the  possibility 
that  the  imperfections  were  symmetric.   Consequently,  a 
least  squares  method  was  used  to  determine  the  best  meridian 
of  data  points  about  which  to  create  symmetry.   This  best 
meridian  was  used  as  the  origin  to   form  the  series 

23   .    . 
W  =  1      Wlnjcos  (n9) 
n=0 

using  the  data  points  from  0  =  0  to  8  =  it  . 
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An  attempt  was  made  to  create  a  finer  meridional  mesh 
of  the  imperfection  data  by  passing  a  general  polynomial 
through  the  fifteen  known  coefficients  for  each  value  of  n 
in  order  to  interpolate  imperfection  coefficients  at  inter- 
mediate stations.   However,  a  polynomial  of  sufficient 
degree  to  pass  through  all  of  the  stations  resulted  in 
unreliable  coefficients  at  intermediate  stations.   Con- 
versely, a  least  squares  fit  using  a  lower  order  polynomial 
resulted  in  a  reasonable  coefficient  representation,  but 
the  coefficients  at  the  original  fifteen  stations  were 
compromised  to  an  extent  that  depended  on  the  degree  of 
the  polynomial.   Similar  problems  were  encountered  using 
a  Fourier  series  representation.   As  a  result  of  the  dif- 
ficulty in  establishing  coefficients  at  intermediate  stations, 
the  fifteen  original  locations  were  selected  as  the  stations 
for  the  numerical  shell  analysis  and  the  conventional  central 
finite  difference  scheme  was  used  to  obtain  the  slope  of 
the  imperfection  shape  away  from  the  boundaries.   At  the 
boundaries  the  simplest  forward  and  backward  differencing 
schemes  were  used. 
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V.   AN  EIGENVALUE  SOLUTION 

To  provide  a  check  on  the  validity  of  the  solution 
from  the  computer  program,  buckling  of  a  simply  supported 
perfectly  circular  cylindrical  shell  was  formulated  as  an 
eigenvalue  problem  and  solved  analytically.   The  set  of 
uncoupled  field  equations  can  be  given  in  the  form 

[E(n)]  {Z(n)"}+  [F(n)]  {Z(n)'}+  [G(n)]  {Z(n)}  =  {e(n)}  (4) 
where 

Z  =  {U(n),  V(n),  W(n),  M^n)}T. 

The  equations  for  the  values  of  the  elements  in  the  [E    ] , 

[F    ],  [G    ],  and  {e    }  matrices  are  given  in  Appendix  C, 

For  n  =  0 ,  the  linear  membrane  solution  N    =  -N   is 

s        o 

assumed.   Thus  the  buckling  modes  of  the  shell  for  the 
boundary  conditions  given  in  equations  (2)  (boundary  condi- 
tion A)  and  (3)  are 

U     =  C-.  cos  iutts/L 

V(n)    =   C2    sin   itnrs/L 

(n)  (5) 

W  =   C-.    sin   mirs/L 

M  =   C.    sin   miTs/L 

s  4 

where   n   ^    0,    C-,    -    C.    are    arbitrary    constants,    m   denotes    the 
axial   mode   number   and   L   is    the    length   of    the    cylinder. 
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Substituting  equations  (5)  into  equations  (4)  results  in  an 
eigenvalue  problem  of  the  form 

[A]  {C}  =  -Nq[B]{C} 

where  [A]  and  [B]  are  4x4  matrices  defined  in  Appendix  C 
and 

{c}  =  {c1/c2,c3,c4}T. 

Because  B  is  singular,  the  inverse  formulation 

~-  [A]{C>  =  -[Bj{C>  (6) 

o 

is  used. 
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VI.   SHELL  DESCRIPTION  AND  INPUT  DATA 

The  shell  used  in  this  study  was  cylinder  A-8  of  Ref.  3 
Its  dimensions  and  properties  are: 

Radius,  r/Re  =  4*°  in« 


Nominal 
length 

Wall 
thickness  ,h 


25  in. 


=  0.00464  in. 


Young's         lc  _    ,A6  , ,  ,. 
modulus, E    =  15.2  x  10   lb/in. 


Poisson1 s 
ratio ,  v 


=  0.333 


Siv  harmonics  of  the  initial  imperfection  data  were  used 
as  input  to  the  program.   They  are  n  =  0,  9,  10,  11,  12,  and 
13.   The  selection  of  these  particular  harmonics  was  based 
on  the  eigenvalue  solution  for  the  simply  supported  shell 
that  predicated  a  minimum  buckling  load  when  m  =  1  and  n  =  9, 
and  the  experimental  observation  in  Ref.  [3]  that  n  =  13 
appeared  to  be  the  critical  mode  where  m  and  n  are  the  mode 
shapes  in  the  meridional  and  circumferential  directions 
respectively.   The  number  of  meridional  stations  used  was 
fifteen  because  no  reliable  interpolation  of  the  initial 
imperfections  could  be  developed,  as  discussed  above.   The 
length  of  the  shell,  L,was  taken  as  8  in.  with  the  first  and 
the  last  data  point  stations  assumed  to  be  the  stations  at 
the  edges  of  the  shell.   Convergence  criterion  was  taken 
as  0.01. 
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As  presented  in  Appendix  A,  the  imperfection  input  data 

_/\       /_  \ 
consists  of  the  rotations  <t>     and  6  '    where 

Sk        9k 

-r(n)  °o   fdW(n)\     _  x(n)    aao  „(n) 

*sv  ='E~     di"   L   and   V   =  ^~       k 
k      o  y     /k         k      o 

in  which  a   is  a  reference  stress  level,  and  E   is  a 
o  o 

reference  elastic  modulus,  a  is  a  reference  length,  and  k 
denotes  the  station  number.   The  data  are  introduced  into 
the  program  through  the  subroutine  IMPERF.   The  FORTRAN 
statements  for  IMPERF  are  given  in  Appendix  D.   The  imper- 
fection data  was  provided  to  the  author  by  the  writers  of 
Ref.  [3]  in  the  form  of  a  deck  of  cards  with  the  data  W/h. 

In  order  to  obtain  a  buckling  load  that  was  a  very 
close  approximation  to  the  biiur cation  lead  or  the  eigen- 
value problem,  the  imperfections  were  reduced  by  a  factor 
10   .   This  technique  has  been  used  successfully  with  this 
program  by  Stillwell  [4]  who  used  very  small  asymmetries 
in  the  load  as  the  triggering  mechanism. 
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VII.   RESULTS  AND  DISCUSSION 

The  solution  to  the  eigenvalue  problem  given  by  equa- 
tion (6)  resulted  in  a  minimum  buckling  load  for  the  simply 
supported  cylinder  of  53.02  lb/in.  when  m  =  1  and  n  =  9. 
The  solution  from  the  computer  program  for  the  minute  imper- 
fections resulted  in  a  buckling  load  of  54.0  lb/in  for 
boundary  condition  A,  simply  supported  edges,  and  displace- 
ment in  harmonics  n  =  0  and  n  =  9.   Some  difference  between 
the  two  results  is  to  be  expected  because  the  analytical 
solution  neglects  the  pre-buckling  deformation  due  to  edge 
constraint,  and  the  numerical  solution  used  only  fifteen 

r-     • .J-.1 

The  computer  results  for  the  simply  supported  cylinder 
with  boundary  conditions  A,  B,  and  C,  with  both  minute  and 
full  imperfections,*  are  presented  in  Table  I.   These  results 
show  a  relationship  between  the  imperfection  sensitivity  of 
the  shell  and  the  boundary  conditions.   For  example,  boundary 
condition  A  shows  a  15%  reduction  in  the  load  carrying  capa- 
city when  the  full  imperfection  is  introduced.   Boundary 
condition  B  shows  a  10%  reduction,  while  condition  C  reduces 
the  capacity  by  4  2.7%. 

The  computer  results  for  the  clamped  shell  are  given  in 
Table  II.   The  experimental  buckling  load  of  shell  A-8  was 


*  -6 

minute  =  full  imperfections  x  10 
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TABLE  I 
SIMPLY  SUPPORTED  CYLINDER 


Loadi: 
Condit 

ng 
ion 

Degree  of 
Imperfection 

Buckling  Load 
(lb/in) 

A 

minute 

54.0 

A 

full 

45.9 

B 

minute 

69.0 

B 

full 

62.7 

C 

minute 

78.6 

C 

full 

45.0 

TABLE  II 

CLAMPED  CYLINDER 

Loading 
Condition 

Degree  of 
Imperfection 

Buckling  Load 
(lb/in) 

A 

minute 

7.5 

A 

full 

no 
convergence 

B 

minute 

43.87 

B 

full 

46.9 

C 

minute 

85.7 

C 

full 

49.1 
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found  to  be  32.87  lb/in  [3].   From  Table  II: 

1)  The  buckling  load  for  condition  A  with  minute  imper- 
fections is  7.5  lb/in.* 

2)  No  converged  solution  could  be  found  for  condition 
A  with  full  imperfections. 

3)  The  buckling  load  for  condition  B  is  slightly  higher 
with  full  imperfections  than  for  the  minute  case. 

4)  The  buckling  load  for  loading  condition  C  is 
decreased  by  42.8%  when  the  full  imperfections  are 
introduced. 

Thus ,  the  stability  of  the  clamped  shell  under  condi- 
tion A  is  drastically  reduced  below  that  of  condition  B  and 
C.   This  is  very  puzzling  and  much  effort  was  expended  in 
search  of  an  explanation.   Unfortunately  -  nnnr  was  round. 
The  computed  result  for  condition  C,  although  generally 
higher  than  the  empirical  result,  shows  the  expected  load 
reduction  due  to  the  initial  imperfections  [3].   Some  of 
the  discrepancy  in  magnitude  could  be  due  to  the  use  of  an 
effective  length  that  was  0.25  in. less  than  the  actual 
length  of  the  shell  since  the  greater  length  would  lower 
the  computed  buckling  load.   The  remainder  may  be  due  to 
the  fact  that  the  actual  boundary  condition  is  some  combin- 
ation of  conditions  A  and  C.   Another  consideration  is  that 
the  symmetric  representation  of  the  actual  imperfections 
may  also  degrade  the  validity  of  the  solution. 


* 

For  this  boundary  condition,  the  imperfections  were 

reduced  by  a  factor  of  10~9. 
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The  circumferential  harmonic  n  =  10  was  found  to  exhibit 
the  most  rapid  growth.   Figure  2  displays  N    vs  W     for 
boundary  condition  C.   The  sudden  decrease  in  slope  occur- 
ring just  prior  to  reaching  N   .    illustrates  the  buckling 

LI  1  L-   • 

(non-convergence)  phenomenon. 
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Nlo)  vs    W(I0)    station     9 


90 
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6C 


70 


N 


0 
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CO 


50 


CLAMPED 

4 
WIO 


40 


W-IO 


CLAMPED 


f* 


30 


10 
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VIII.   CONCLUSIONS 

1.  The  imperfection  sensitivity  of  a  cylinder  is  dependent 
on  the  boundary  conditions . 

2.  The  computer  program  yields  more  realistic  results  for 
the  simply  supported  boundary  condition  than  for  the 
clamped  condition. 

3.  The  initial  imperfections  result  in  a  significant 
reduction  in  load  carrying  capacity  in  the  simply 
supported  case  but  less  reduction  than  expected  [3]. 

Further  investigation  is  necessary  to  determine: 


experimental  results? 

2.  Why  is  the  clamped  boundary  condition  A  so  sensitive  to 
imperfections  in  the  computer  program? 

3.  Why  does  loading  condition  R,  clamped,  appear  to  be 
strengthened  by  the  introduction  of  full  imperfections 
in  the  computer  solution? 
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APPENDIX  A 

This  section  presents  the  appropriate  equations  for 
including  initial  imperfections  in  the  computer  program. 
The  strain-displacement  relationships  used  in  this  analysis 
are 

e    =  U1  +  W/R  +  ($2  +  $2)/2 
s  s      s 

ee      =   V/r   +    r'U/R  +   W/R     +    ($2  +  $2)/2  (Al) 

c        =    (V  +  u'/r    -    r'V/r   +    $    $fi)/2 
where 

<*       =  -V71   J-   TT  /ID 

s  s 

$e   =-W/r  +  V/R0  (A2) 

9    =  (V  +  r'V/r  -  U*/r)/2 
and  where  e  ,  eQ,  e  .  represent  Fourier  coefficients  of  the 

S     0     StJ 

reference  surface  strains,   $ ,  C>  ,  and  $  „  represent  refer- 

6        s6    c 

ence  surface  rotations,  and  Rfi  and  R   are  the  principal 
radii  of  curvature.   The  superscript  primes  and  dots  denote 
partial  differentiation  with  respect  to  s  and  G  respectively 

Consider  an  imperfect  cylinder.   Before  loading,  the 
imperfections  appear  as  a  displacement  from  the  perfect 
cylinder.   Note  that  this  is  a  displacement  with  no  accom- 
panying strain.   Call  this  displacement  W,  and  the  total 
displacement  from  the  perfect  shape  W* .   Then 
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w*  =  w  +  w 

where  W  is  that  part  of  the  total  displacement  caused  by 

the  applied  load  and  is  associated  with  strain.   Replacing 

W  with  W*  in  equations  (Al)  and  (A2)  and  requiring  that 

£   =£n  =  e,.  =  0  when  U  =  V  =  W  =  0,  leads  to 
s     6     s6 

E   =  U'+W/R  +  [(-W'+U/R  )2  -  2W1  (-W+U/R  )  +  $2]/2 

eQ  =  V/r+r'U/r  +  W/RQ+  [(-W'/r+V/R  )  2  -  2W"  /r  (-W  /r+V/Rg) 
+  $2]/2 


e   =  [V'+U'/r  -r'V/r  +  (-W+U/R  )  (-W/r+V/Rfi)  -  W 
sy  s  w 

(-W/r+V/Rfi)  -  W'/r(-W  +U/R  )]/2 


Defining  the  imperfection  rotations  *   and  <!>  »c 


$   =  -W1      and      $Q  =  -W/r 

s  9       ' 


equations    (A2)    become 


e      =   U'+W/R      +    (3>2+2$    $     +$2)/2 
s  s  s  s    s 

E      =   v*/r  4-r'U/r   +   W/Ra  +    ($2  +2?„$fl+  $2)/2 


e      =    (V'+U'A   -    r'V/r   +    $    $.  +  $    $.  +  $    $Q)/2 

se  sesese 


From  Ref.     [2] 


-°  Z      <f>(n)cos    n6    2    =   =2      E      6(n)cos    n( 

Eol      n=0      S  Eo   n=0      S 


Redefine    3      such    that 
s 
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=£     Z      3(n)    cos   n0   =   $2   +    2$   $ 
Eo  n=0      s 


.(n)„    n6]2  +  2f 


=  E"   \\    ZA      cos 
o      \n=0 


\n=0 


Z      <J>        cos    n6         Z      if)        cos    ni 


n=0 


Similarly 


^      Z      B6    cos    n0    =    $q+  2$ 
o  n=0 


^  <        Z      (Pgcos    ne]2  - 


00  \    /  00 

Z      ^D(n)sin   n9         Z      (fv(n) 


n=l 


n=l 


and 


e2-        En   6sesin  ne  =  V 

o        n=l 


+    $$+$< 
s    6  s 


f  O     \  o    I  /     °°  .  1 1      °° 


in  i 


sin   no      + 


Z      <t>         cos    n( 
n=0      S 


Z      <J>  2      sin   n( 

T  O 

n=l 


Z      (J)  (n)cos   ni 
n=0      S 


Z      <}>        sin  nl 
n=l       L 


where 


Ys 


dw 


(n) 


d£ 


and 


r(nJ         nw(n) 


(n) 


.   P  ,  t  a   -(n)    =  E   Wv    7/(aa    )  . 

and   t,    =    s/a,    p    =   r/a,    and   w  o  o 


The    changes    to    the   program  are    in   SUBROUTINE   PHIBET 
given  below. 
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SUB 

COM 

1/IB 

2/IB 

3/IB 

3/IB 

4,10 

6/IB 

2/IB 

8/BL 

O/BL 

9/3L 

3/BL 

1/BL 

2/BL 

3W3( 

4/Bl 

COM 

COM 

OX  = 

OT= 

RRA 

GA= 

KP2 

DO 

EN  = 

IK= 

U3( 

V3( 

W3( 

PHI 

PHI 

1  PHI 

IF( 

KP= 

PHX 

60  PHT 


ROUTINE 
HON 

L1/MNMA 
L2/N( 10 
L4/KMAX 
L7/MNMA 
)  ,  J  0  (  1 0 
L 12 /KM A 
L13/ITP 
6/Z(4,2 
8/R(20C 
10/ PHI  X 
1  I/O* XI 
12/TOLI 
15/NU,U 
10) 

27/BX3( 
MON/SLI 
MON/BLP 
CMXI (K) 
CMT(K) 
=1./R(K 
GAY(K) 
=  K+2 
1  M=1,K 
N  {  M ) 
KP2+(M- 
M)=Z( 1, 
M)=2(2  i 
M)=Z (3, 
X(M)=-T 
T(M)=EM 
(  H )  =  ( 1  D 
ITRMAX. 
DO  6 
K+(f'.-l  ) 
IM»=PHI 
(M)=PH1 


PHIBET (K) 


X 

) ,MNINIT 

»KL 

XO,M 

,10) 

X 1 ,  K 

MAX., 

2  CM  , 


AXD(  1 
,  IJS( 


10  ) 

200 
TOP 


N-AX2, 
LS?-AX 
SOE,0 
M(20C 
,PH  IT 
)  ,PHE 
L 


0),MAXS(1G) ,MAXSY(10),IS(10,10),JS(1G,10),ID(10 

10) 

NCONV 

SE, ALOAD 
) ,0MT{200) 

(10) t PHI ( 10) 
E,T0,T2 


1(10) ,V1 (10) ,W1(10) ,V2(1C) ,U2(10) ,W2(10) ,U3( 10) ,V3(1C) , 

10) ,BT3< 10) , BXT3(1C) ,BE3< 10) 
MP/PHIXB<200) , PHITB(2C0) 
HS/PHX( 10) ,PHT(1C) 

) 


1) 
IK 

IK 
IK 
DL 
*W 
LI 

EC 
rj 

*K 
XP1 
TB 


*KMAX 

) 

) 

) 

I*(W3 

2  (  M )  * 

*(V3< 

•  1  )  R 

N=l  ,M 

MAX 

(KP) 

(KP) 


(M)-wi(M)  )+0X*U2(M) 
RRA+V2<M)*0T 

;')-Vl  (M)  )+GA*V2(M)+EN*U2<  M)*RRA)*.5 
ETURN 

NMAXO 


pn  q  v=i 


I  m  a  y 


+PHX( J)*DHIX(I ) 


SMT=0. 

s:-'.R=o, 

SMF=0. 

IF(N('-' ).EO.C  )  GO  TO  20 

MAXL  =  iv.AXS<  Ml 

IE( KAXL.EO.O)  GO  TO  2 

DO  ?•  L=l  ,  MAXL 

I  =  IS(L,",) 

J-JS(LtM) 

SMO=SMO+PHIX(I J*PHIX( J) 
1  +  PHXU  )*  PH1X<J) 

SMT=SMT-PHIT( I  »*PHIT( J  J 
1  -  PHTd  )--:-PHIT(j  )  -PHT{ J)*PHIT( I ) 

SMR=SMR+PHIX(I )*PHIT( JJ+PHIX ( J )*PHIT( I ) 
1    +  PHX(  T  }*PHIT<  J  )+PHX(  J)*PHIT<  I  )+PHIX(  I  )*PHT(  J)  +  PH  I  X  (  J  )  *PHT  (  I  ) 
3  SMF  =  SMF-P'ri!  (  I  )*PH1  (  J) 
2  MAXIMA/DP'*) 

IF(MAXL.EC.O)  GO  TO  4 

DO  5  L=1,MAXL 

I=ID(L,r) 

j=jp(l,:')  .  . 

SM0  =  Si10+PHIX(I  )*PHIX( J  ) 
1  +  PHX( I )*PHIX( J) 

SMT=SMT+PHTT( I  )-PHIT( J) 
1    +    PHT(I)*PHIT<J) 


+PHX( J) *PKTX(I ) 


+  PHT(  J)--'PHIT(  I  ) 
SMR=SMP-PHIX<  I  )*PHIT(  J  )+PHIX(  J)--PHIT(  I  ) 
1       -PHX( I ) vPHIT (J )+PHX( J)*PHIT( I )-PHIX( I )*PHT( J)+PHIX(J)*PHT(I) 
5    SMF=SMF+PHI (I)*PHI ( J) 
4    IF<KAXSY(M).EQ.C)     GO    TO    10 
1=1 JS(M) 
SMC=SMC+PHIX( I )**2/2. 
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!0 


!1 


1 


LG 


+  PHX( I )*PHIX(I  ) 
SMT=SMT-PHIT(I )**2/2. 

-PHT( I)*PHIT(I  I 
'SMR=(SMR+PH1X( I)*PHIT( I) ) 

+PHX( I )*PHIT< I )+PHIX( I ) 
SMF=SMF-PHI  (I  }**2/2. 
GO  TO  1C 

DC  21  L=ltMNMAXO 
SMC=SMC+PHIX(L)**2 

+2.VPHX (LJ-PHIX(L) 
SMT=SMT+PHIT(U**2 

+2.*PHT(L)*PHIT(L) 
SMF=SMF+PHI  (L  )  -'--2 
IFCM.GT.MMMAXOl  GO  TO  11 
SMG=SM0+PH1XIM)**2 

+  2.*PHX(f-'.}*PHIX(M) 
BX3(N)=SKO*.5 
BT3(H)-SMT*.5 
BE3(M)-SMF*.5 
BXT3(M)=0. 
GO  TO  9 
BX3( M)=SHO 
BT3(f';)=SMT 

BXT3(M)=StfR*.5 
BE3(K)-SMF 

CONTINUE 

RETURN 
END 


PHT(  I) 
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APPENDIX  B 

The  load  applied  in  the  computer  program  is  uniform  and 
taken  in  the  first  mode  only,  i.e.,  n  =  1.   To  accomplish 
this,  the  following  changes  were  made  to  the  program  logic: 


SUBROUTINE  ZANDZ  531 


DO  15  J  =  1,  4  659 

*    IF(M.NE.l)  ELLS (J)  =  0 


RETURN 
END 


SUBROUTINE  FORCE (K)  960 


DO  21  J  =  1,  4  081 

*    IF(M.NE.l)  ELLS (J)  =  0 


RETURN 
END 


*  = 


Required  inserts 
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APPENDIX  C 

The  elements  of  the  [A]  and  [B]  matrices  used  in  the 
eigenvalue  solution  are  defined  below. 

All  =  Gll  "  E11q2  A21  =  -F21Q 

A12  =  F12Q  A22  =  -E22Q2  +  G22 

AX3  =  F13Q  A23  =  G23  +  E23Q2 

A14  =  °  A24  -  G24 


A31  =  -F31Q  A41  =  ° 

A32  "  G32  "  E32°2  A42  =  G42 

A,,  =  G,,  -  E,,Q2  A,,  =  G,,  -  E„,Q2 


33  33     33"  43     43     43 

2 

34  =  G34  "  E34Q  A44  ~  G44 


B13  =  1/4  n2  B21  =  1/4  nQ 

B12  =  1/4  nQ  B22  =  1/4  Q2 

B13  "  °  B23  "  ° 

B14  =0  B24  =  0 


B31  =0  B41  -  0 

B32  =0  B42  =  0 

B33  =  Q2  B43  =  0 

B34  =  °  B44  "  ° 
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where   Q  =  niTs/L 


and 


'11 


'12 


J13 


'14 


b 

0 
0 
0 


'21 


=  0 


E 


E 


b(l-v)    d(l-v)  |  3  \2 


22 


23 


(1-v)  3n 
2r    RA 


'24 


=  0 


Re 


E31  =  ° 


E 


&)iu 


'32 


'33 


'34 


11 


12 


13 


14 


=  E 


23 
d(l-v) 


2n 


=  1 


=  0 


2r 


8r 


=  bv-d 


=  0 


(l-v)n' 
2r2Rn 


41 


'42 


'43 


'44 


0 
0 

-d 
0 


(l+v)bn  +dn(l-v)    -3 


R, 


21 


22 


23 


24 


=  -F 
=  0 
=  0 
=  0 


12 


31 


32 


33 


34 


=  -F 
=  0 
=  0 
=  0 


13 


41 


42 


43 


44 


34 


and 


'11 


(l-v)bn' 

o  2 
2r 


+  d(l-v) 


-1 
8Rer: 


n 


G12  =  0 


G13  =  ° 


G14  =  ° 


21 


=  0 


'22 


=  0 


23    rR 


+  (dn(l-v) 


2r 


-2  (l+v)n' 
r2R, 


-vn 


'24 


31 


'32 


'33 


'34 


'41 


'42 


'4  3 


=  0 


Z|n  +  dC^vin  (2(1+v)| 
-b    d(l-v)n2 


r2R, 


(1+v) 


■n- 


Re 


vn 


0 

dvn 

rRe 

dvn2 


'44 


=  1 
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APPENDIX  D 

The  following  subroutine  was  used  to  calculate  <j) 

ana  *<">  in  tne  colter  prolan,. 
k 


SUBROUTINE  IMPERF 
DIMENSION  FNT(15) ,DFNT(15) 
COMMON  /BLIMP/PHITB(200) ,PHIXB(200) 
A0=4 

SIGA=1000. 
EO=15200000. 
H=8./14. 
Hl=l./H 
H2=1./(2.*H) 
DO  1  1=1,6 
AJ2=I+7 

IF(I.NE.l)  GO  TO  8 
AJ2=0. 
8  REAP ( 5 ,  LO )  ( FNT ( J )  . .  I  =  3  . ] b  ) 
DO  2  K=2,14 

2  DFNT(K)=(FNT(K+1)-FNT(K-1) ) *K2 
DFNT(1)= (FNT (2) -FNT (1) *H1 
DFNT(15)= (FNT (15) -FNT (14) ) *H1 
KR=I-1 

DO  3  Kl=l,15 

KZ=KR*15+K1 

PHITB (KZ)  =  (FNT (Kl) *AJ2*EO/SIGA) *  (0. 00464) 

PHIXB(KZ)=(DFNT(K1) *EO/SIGA) * (-0.0  0  464) 

3  CONTINUE 
1  CONTINUE 

10  FORMAT (6E12. 5) 
RETURN 
END 
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boundary  conditions  of  the  cylinder. 
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